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The transition from a flowing to a static state in a granular material is studied using large-scale, 
3D particle simulations. Similar to glasses, this transition is manifested in the development of a 
plateau in the contact normal force distribution P(f) at small forces, along with the splitting of the 
second peak in the pair correlation function g(r), suggesting compaction and local ordering. The 
mechanical state changes from one dominated by plastic intergrain contacts in the flowing state to 
one dominated by elastic contacts in the static state. We define a staticity index that determines 
how close the system is to an isostatic state, and show that for our systems, the static state is not 
isostatic. 



Granular materials are many-body systems which can 
exhibit properties of both liquids and solids, often at the 
same time. Features such as jamming in silos and hop- 
pers |lj and heterogeneous force propagation || have 
motivated new proposals for constitutive relations for 
granular materials \pw%, in contrast to the established 
elasto-plastic theories Mp| . Moreover, the emerging field 
of "jammed systems" |7| seeks to understand whether 
analogies may be drawn between diverse systems which 
have the common properties that they are far from equi- 
librium, and unable to explore phase space. 

O'Hern et al. || have recently shown that similarities 
exist between the properties of static granular packings 
and other amorphous systems; in particular, they com- 
pared force distributions for simulated systems of soft- 
spheres undergoing a glass transition with experimental 
data from static granular packings ||||. These studies 
suggest that it may be interesting to consider the tran- 
sition of granular materials from flowing to static, and 
inquire about its similarities with the glass transition. 
In this Letter we study the dynamic jamming transition 
of systems of athermal grains through large-scale simu- 
lation. In particular, we consider the behavior of dense 
packings of granular particles flowing down an inclined 
plane as the tilt angle 8 is reduced through the angle of 
repose 6 r at which flow ceases. 

The chute flow geometry we employ relies on gravity 
to compactify the system. Our earlier simulation stud- 
ies in this geometry [ fU3|]IT| have proven their reliability 
by reproducing key experimentally observed characteris- 
tics of dense granular flows [Q , such as the existence of 
steady state flow over a range of angles between r and a 
maximum angle of flow stability # max , as well as the de- 
pendence of the angle of repose on pile height for shallow 
piles less than 20-30 grain diameters in depth. 

The transition from flowing to static states suggests 
that for systems generated in this way, 6 r can be regarded 
as an analogue for the glass transition temperature T g , 
and we characterise the transition by the following prop- 
erties: 



• P(f), the probability distribution of normal forces, 
develops a plateau in the static state. 

• The radial distribution function, g(r) develops a 
split second peak, and growth of the first peak, in- 
dicating compaction and increased short-range or- 
dering. 

• The average co-ordination number, z c , jumps dis- 
continuously. Moreover, the static state that re- 
sults when motion ceases is not isostatic. 

Furthermore, we observe a significant change in the 
number and nature of interparticle contacts: although a 
finite fraction of contacts are plastic (at Coulomb yield) 
for all flowing piles, almost none are plastic when the pile 
becomes static The average co-ordination number 

increases continuously towards four in the flowing state 
and jumps to well above four for the static packing - four 
being the theoretical minimal co-ordination number for 
a network of frictional grains. 

We carried out molecular dynamics simulations in 3D 
on a model system of N = 8000 mono-disperse, cohesion- 
less, frictional spheres of diameter d and mass m. The 
system is spatially periodic in the xy (flow - vorticity)- 
plane, and is constrained by a rough bottom bed in the 
z-direction, with a free top surface. The static packing 
height is about AQd. Particles interact only on contact 
through a Hertzian interaction law in the normal and tan- 
gential directions to their lines of centers |l4||L5||. That 
is, contacting particles % and j positioned at and in- 
experience a relative normal compression 



where rjj = — Vj , which results in a force 

F y =F„+F t . (1) 
The normal and tangential contact forces are given by 

/ Tfi \ 

F„ = V S/d ( k n Sxiij - — 7„v„ I , (2) 
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Ft - (-ktAst - y7tvt) , (3) 

where ny = ry/ry, with ry = |ry|, v n and v t are the 
normal and tangential components of the relative surface 
velocity, and fc„ jt and 7 n> t are elastic and viscoelastic con- 
stants respectively. As t is the elastic tangential displace- 
ment between spheres, obtained by integrating surface 
relative velocities during elastic deformation of the con- 
tact The magnitude of As t is truncated as neces- 
sary to satisfy a local Coulomb yield criterion F t < [iF n , 
where F t = |F t | and F n = |F„|, and fi is the local particle 
friction coefficient. 

In the steady state flow regime, changes in the par- 
ticle friction /x, the damping coefficients j n t, or parti- 
cle hardness k n , do not change the qualitative nature 
of our findings jllj. For the present simulations we set 
k n = 2 - I0 5 mg/d, k t = ~k n , J n = 50y^gjd, j t = 0.0, and 
H = 0. 5 |p7|| . We choose a time-step St = 10~ 4 t, where 
T = y/djg. 




FIG. 1. Bulk packing fraction cf> and bulk averaged kinetic 
energy (KE) per particle, normalised by mgd, as a function of 
tilt angle 8. The dashed arrow indicates the value of <f> when 
6 is reduced from 19° to zero. 

To initiate flow, we tilt the bed to a large angle 
[6 « 26°) long enough to remove any history effects of 
construction. We then lower the angle, in unit incre- 
ments, to the desired value in the range for which there 
is steady state flow, 20° < 9 < 26° |uj. In this range, the 
energy input from gravity is balanced by that dissipated 
in collisions. We can approach the flowing angles from 
above and below and always obtain the same results. 

Lowering 9, the kinetic energy of the system decreases, 
however flow continues down to Q r {— 19.4° for the set of 
parameters used here Jl^]). As 9 is reduced, the vol- 
ume fraction increases, as shown in Fig. [|. We generated 
further static packings by taking the stopped state at 
9 = 19° and reducing the tilt to zero. The volume frac- 
tion of the static packing for 9 = 0.0, = 0.599 ± 0.005, 
is less than the random close packing value 4> rC p — 0.64. 



In order to investigate the distribution of contact forces 
in flowing and static states, we computed the probability 
density function (PDF) P(f), where / = F n fF n t z \ is 
the ratio of the normal contact force magnitude F n to 
its average value at the depth z of each contact. Since 
our system is spatially periodic in the horizontal plane, 
with no side- walls, the pressure does not saturate with 
depth, requiring separate averaging at each value of z. 
The data were averaged over 1000 configurations over a 
period of 200t in the steady state at each flowing angle 
and 5 configurations for the static states. 

As shown in Fig. the PDFs exhibit familiar exponen- 
tial tails at high forces. When 9 is reduced towards 9 r , 
a plateau develops near / = 1 (see inset), similar to the 
behaviour observed in Ref. near the glass transition. 
It may be argued on this basis that 9 r appears analogous 
to T g in a glass-former. However, the development of 
this feature is rather subtle in our system, and does not 
appear to be a sharp indicator of the transition from the 
flowing to the static case. Upon further relaxation of the 
stress by reducing 9 to zero, the feature becomes more 
prominent similar to the results in Ref. H. 




FIG. 2. Distribution of normal contact forces P(f) for var- 
ious tilt angles (see text). The inset shows the emergence of 
a plateau in P(f) at small forces for the static systems. 

As with glass transitions, we have found only subtle 
structural evidence for the transition from flow to rest. 
The radial distribution function g(r) is fairly insensitive 
to the tilt angle, as shown in Fig. ||. Although there is 
no long range ordering, in the inset to Fig. ^ we observe 
the gradual splitting of the second peak as 9 is reduced. 
The system develops more locally ordered structures as 
it becomes more compact |l9| . 

More significantly, the first peak in g(r) increases with 
decreasing 9. Because particle neighbors are defined only 
on contact in our simulations, the first peak in the ra- 
dial distribution function corresponds to the average co- 
ordination number z c , which correspondingly increases 
as 9 decreases (see Fig. ||). 
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FIG. 3. Radial distribution function g(r) for various tilt 
angles 9. The inset shows the region near the second peak. 



An important characteristic that separates the granu- 
lar system from a glass forming liquid is the existence of 
interparticle friction and yield. To explore the influence 
of intergrain friction, we distinguish between intergrain 
contacts at their yield criterion F t = fiF n (plastic con- 
tacts), which respond plastically to some external pertur- 
bations, from contacts with F t < \xF n (elastic contacts), 
whose response is determined by elastic deformations of 
the constituent grains. Note that the relationship be- 
tween plasticity/elasticity of individual contacts between 
grains and the mechanical response of the entire system 
has not been explored in this study. 

We find that the nature of grain-grain contacts is quite 
different in the flowing and static regimes. The proba- 
bility densities of frictional saturations £ = F t / \xF n are 
shown in Fig. ^ for various tilt angles. The fraction n c 
of plastic contacts (with £ = 1) decreases as the tilt an- 
gle is decreased towards 9 r (see Fig. ||) , accompanied by 
a decrease in the average frictional saturation of elastic 
contacts. For the static piles at 9 < 6 r , almost all grain 
contacts are elastic. 

A related question is whether or not static piles of rigid 
grains satisfy an isostaticity condition, where the num- 
ber of contacts is the minimum required to satisfy force 
and torque balance equations for each grain po|| . This 
isostaticity hypothesis has been frequently invoked in re- 
cent theoretical work attempting to derive macroscopic 
constitutive relations directly from the microstates of 
static granular packings plj . Isostaticity is required for 
a unique determination of the individual contact forces 
solely in terms of the microstate. Otherwise, the details 
of grain deformations (for example, Hookean or Hertzian 
force laws) must be considered in order to determine the 
stress state of the pile p3|, as the rigid grain problem 
becomes ill-posed. 




FIG. 4. As the tilt angle is decreased towards 9 r , the 
fraction of plastic contacts n c (squares) decreases while the 
co-ordination number z c (diamonds) and the staticity index 
s c (circles) increase. The isostatic values of z c and s c are 
indicated by the solid arrows. The dashed arrows indicate 
the values of z c and s c when we take 9 down to zero from 
a stopped state. The point indicated by * at 9 = 20° is for 
k n =2- 10 7 mg/d. 



Counting the degrees of freedom (DOF) in the con- 
tact forces suggests that isostaticity requires z c = 4 for 
frictional contacts (with 3 DOF - one vector force - per 
contact and 6 equations per grain) and z c = 6 for friction- 
less contacts (with 1 DOF - one normal force - per con- 
tact but only 3 equations per grain) for spherical grains. 
However, one DOF is eliminated for each frictional con- 
tact that reaches the yield criterion and becomes plastic 
p2fl : the magnitude of the tangential force in this case 
is determined by the normal force and the friction coef- 
ficient [i. Thus, if the fraction of plastic contacts is n c , 
for a co-ordination number z c the total number of DOF 
characterizing interparticle forces is 



N 

N f = -J Zc t^ 1 ~ ™ c ) + 2n d 



(4) 



and requiring Nf > 6N to match the number of kine- 
matic constraints on the particles yields 



(3 - n c )z c /2 > 6, 



(5) 



where we have defined the staticity index s c . The equal- 
ity corresponds to isostaticity. By analogy to the sys- 
tem in Ref. |22j , the response of "hyperstatic" piles with 
s c > 6 to external forces cannot be determined without 
considering the history of grain contacts. 

Figure [| depicts s c as a function of 8. As expected, 
s c < 6 for all flowing piles, though it increases with de- 
creasing 9. However, it appears that a finite jump of s c 
at 9 r leaves the static pile significantly hyperstatic. Note 
however that the point denoted by * in Fig. [| is for a 
system where k n — 2 ■ 10 7 mg/d, indicating that grain 
hardness is a relevant parameter in determining s c and 
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FIG. 5. The probability distribution P(C) of frictional satu- 
ration £ = Ft I fiF n at various tilt angles. In the flowing state 
a finite fraction of plastic contacts exist (see Fig. ^). 



What do these observations teach us about the nature 
of the "jamming" transition of the granular medium? 
Our earlier work showed that granular flows obey a local 
rheology independent of the history of loading M,0|. 
Also, static piles are known to generally be able to sup- 
port a range of external stresses without any macroscopic 
rearrangement. So, how was the original stress state and 
associated configuration of grains in this static pile deter- 
mined? One possible answer is that this state was deter- 
mined by the local loading conditions at the instant the 
flow ceased. This resolution is akin to ideas discussed in 
the fixed principal axes (FPA) model m , except that the 
stress state of the static pile, including principal stress 
axes, do change subsequent to flow arrest (which is man- 
ifested as burial in the FPA model) as a function of ex- 
ternal loading, due to contact elasticity of grains. 

We have observed analogies between jamming of gran- 
ular media and the glass transition. A plateau in the 
force distributions appears at small forces once the sys- 
tem ceases to flow, as observed in numerical studies of 
liquids at the glass transition M. Also, more locally or- 
dered structures appear as the system compacts. How- 
ever, the most significant effect is the transition from a 
flowing plastic system with a local rheology to a static 
system that may preserve the fingerprint of the stress 
state at the instant of jamming, which makes this tran- 
sition similar to a glass transition, in that glasses also 
have an ability to be pre-stressed. Finally, we note that 
we have performed a similar study with Hookean spheres, 
with essentially the same results. 
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